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We consider a class of deterministic susceptible-infectious- removed 
(SIR) models on time-independent weighted contact networks where 
the infection and removal processes are Poisson. We prove that a 
particular closed form of these models generates an exact representa- 
tion of the expected infectious time series for networks with no closed 
loops. This is the first exact deterministic representation of a complex 
heterogeneous and finite epidemic system which can be numerically 
evaluated. 



1. Introduction. 

1.1. Background. The majority of epidemic models fall either into the 
category of stochastic simulation [2, 4] or into the category of determinis- 
tic differential equation-based models [1, 11]. These two strands developed 
largely independently for much of the twentieth century. Thus, an interesting 
question arises as to the precise mathematical connection between stochastic 
and deterministic models. 

The first epidemic models were based on the assumption that populations 
are evenly mixed, with each individual equally likely to interact with any 
other individual at any time [6]. A classic example of this type of model 
is the Susceptible-Infectious-Removed (SIR) compartmental model whereby 
individuals are classified according to being in one of these three states. It 
has been shown that for this type of mean-field model, the average of many 
stochastic simulations (the expected outcome) converges to the solution of 
the "equivalent" mean-field deterministic model in the limit of an infinite 
population size and subject to strict conditions regarding the initialisation 
of the epidemic [20]. In other words, the expected infectious time series in 
large populations is given by the deterministic representation. 
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More recently, a higher degree of reahsm has been introduced by con- 
sidering stochastic simulation models on contact networks where individ- 
uals are only able to contact a limited subset of the population. This en- 
ables significant heterogeneity to be incorporated, treating individuals as 
distinct entities with fixed connectivity to pre-allocated neighbours. While 
stochastic simulations are readily extended to incorporate such systems, de- 
terministic descriptions have been more problematic. Several methodologies 
have been developed including pair-approximation models [8, 9, 16], degree- 
based models [15], and models based on the probability generating function 
(PGF) formalism which are applicable to configuration networks [22] as well 
as the related edge-based compartmental modelling [13, 14]. It has been ob- 
served [7] that these models are, at some level, equivalent and are all derived 
from similar principles of independence. However, although comparison with 
stochastic simulation can sometimes be good, the link between these models 
and their stochastic counterparts remains obscure. 

There are two idealised scenarios in which exact correspondence between 
stochastic simulation and deterministic descriptions has been shown. Firstly, 
correspondence sometimes occurs in the limit of infinite populations for par- 
ticular idealised graphs [3, 5] which cannot be exactly realised in practice. It 
can also occur with some very simplified systems whose symmetry properties 
can be exploited to achieve reductions in the stochastic description [10, 19]. 

A modelling framework which emphasises the link with stochastic simu- 
lation is the subsystem approach [17, 18]. This splits the epidemic system 
up into smaller components governed by their own exact subsystem pro- 
cess equations. These can then be used to construct deterministic models 
via suitable closure assumptions. Here we consider the example of an SIR 
model on finite networks governed by transmission and removal processes 
which are Poisson. We show that by considering subsystems at the level of 
pairs, a closure can be proposed which determines the expected infectious 
time series exactly for arbitrary network structures with no closed loops and, 
in some special circumstances, for particular networks with closed loops. 

1.2. Statement of the main result. We consider an SIR compartmental 
model composed of P individuals whose states are described at any given 
point in time by vectors I and S with respective components Ij and Si, 
i € {1,2, P} such that /j = 1 if individual i is infectious (Ij = otherwise) 
and = 1 if individual i is susceptible {Si = otherwise). Transmission 
and recovery occur by Poisson processes with rate parameters Aj and i^i 
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respectively where: 

p 

Aj — ^ ^ Tij Ij Si, Hi — Ti-^i) 

where T is a "transmission" matrix with (time-independent) elements Tij 
denoting the rate parameter for an infectious node j infecting a susceptible 
node i [Ta = for all i) and where 7^ denotes the rate parameter for an 
infectious individual i to recover. 

As shown by Sharkey (2011) in [18], for any transmission matrix T and 
any nodes i, j the following differential equations hold: 

P 

i=i 
p 

p p 
{Silj) = ^ Tjk{SiSjIk) - ^ TikihSilj), 

k=l,kj^i k=l,k^j 

-Tij{saj)-j,{sj,), 
p p 

(1.1) (SiSj) = — ^ TikihSiSj) — ^ Tjk{SiSjIk), 

k=l,k=ij k=l,k=/=i 

where {Si) and (/j) denote the time-dependent probabilities for individual 
i to be susceptible and infectious respectively, and expressions of the form 
(AiBj) denote the time-dependent probability that individual i in state A 
and individual j is in state B with a similar interpretation of terms of the 
form {AiBjCk)- Here and throughout, we adopt the dot notation to denote 
time derivatives. 

Note that it is a short step (see Sharkey (2008) [17]) from (1.1) to the 
familiar population-level pair equations [8, 9, 16], also proved independently 
by Taylor et al. (2012) in [21] for the susceptible-infectious-susceptible vari- 
ant. 

This system can be completed by formulating differential equations for 
the triples, quadruples, and so forth until we reach the full system size. This 
yields a self-contained system of differential equations that exactly deter- 
mines the probabilities of each quantity given initial conditions. However, 
cascading these equations up to the full system size will usually result in 
a system that is impractical to solve due to its sheer size. This is why this 
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system is typically closed at some level by introducing a functional relation 
defining higher-order probabilities in terms of lower-order ones. One of the 
most frequently used closure relations can be written as 

(1.2) {A.B,C,) . i^ilM^ 

for the current context. Applying this closure relation to our system at the 
level of pairs we arrive at the following system: 

p 

{X,) = -Y,T,,{X,Y,), 

i=i 
p 

{¥,) = ^r,,(x,y,)-7,(y,), 

_ >f T iWj){X,Yk) ^ {X,Yk){X,Y,) 

k=l,k^i k=l,k=ij 

-Tij{X,Yj)-^j{XiYj), 

>f ^ {YkX,){XiX,) ^ {Xaj){X,Y,) 



(^^) k=i,k^i ^^j) 



k=l,k=/=j \ " k=l,kj^i 



(1.3) 



where we use X for susceptible and Y for infectious to emphasise that these 
are approximating differential equations based on the closure. 

Our main aim is to show that when matrix T contains no closed loops 
(i.e. is tree-like) and the system is in initiated in a definite state (that is, one 
of the 3^ possible configurations has probability 1 at time t = 0), then the 
system can be closed at the level of pairs such that the closure holds exactly. 
Specifically, solving the closed system above, we obtain the same values for 
all probabilities present in the unclosed system: {Xi) = {Si), (Yi) = (/j), 
with similar equalities holding for the pairs. 

In fact, we will prove the following theorem. 

Theorem 1.1. Let us assume the following: 

• The graph (transmission network) is a tree (no loops). 

• The initial condition is not mixed, i.e. the system is initially in one of 
its 3^ possible configurations with probability 1. 

Then the following relations hold: 

{Sj){SiS,h) = {S,Sj){SjSk) 
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for all i G {1,2,..., P} and for all j with links towards i and all k with links 
towards j, i ^ k; 

{s.i){ikSiij) = {ikSi){sa,) 

for all i € {1, 2, P} and for all j and k with links towards i, j ^ k. 

The theorem wiU be formulated in a more general context stating that 
even higher order closure relations are also exact. 

Figure 1 shows the numerical solution of (1.3) for a small network of 9 




Fig 1. a) A simple undirected tree-like network indicating two nodes which we infect to 
initiate epidemics, with all others nodes initially susceptible, b) The mean (dots) of 100,000 
stochastic simulations on the network with transmission rate r = 0.1 across each link and 
removal rate 7 = 0.05 for each node, with errorbars denoting the 5th and 95th percentiles 
plotted together with the solution of (1.3) (solid line) using the Matlab code published with 
[18]. 

nodes where it is clear that it is accurate to within the precision visible on 
the graph. Matlab code for solving the system of equations (1.3) is provided 
in Sharkey (2011) [18]. This code also works on networks with closed loops, 
but is no longer exact in these cases. Closed loops of order three utilise the 
alternative closure {AiBjCk) = {AiBj){BjCk){AiCk) / {Ai){Bj){Ck) which is 
believed to gain increased accuracy in most circumstances, but these do not 
occur in the tree-like graphs considered in the present work. 
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The structure of the paper is as fohows. Section 2 introduces some no- 
tation which is needed to prove the result. This also contains an important 
theorem (Theorem 2.1) which specifies equations describing the probabili- 
ties of the states of arbitrary subsystems (the proof of this result is given in 
Appendix A). The relevant state space for our domain of a tree- like (unclus- 
tered) graph is then developed. Section 3 proves the main result, initially 
focusing on some special cases to help motivate and facilitate understanding 
of the main ideas and steps of the general proof in Section 3.5. The main 
ingredient for the general proof is Lemma 3.1 which is proved via Theo- 
rem 2.1. Theorem 3.1 then follows easily by induction from Lemma 3.1. The 
theorem as stated above is a simple corollary of Theorem 3.1. In Section 4 
we discuss an application of the pair-based model to some special cases of 
graphs with loops where it is also exact. 

2. Formulating the full system. In this section we introduce a new 
notation which will assist in formulating the set of differential equations for 
the full system. In (1.1) we formulated the diff'erential equations up to the 
level of pairs and we said that this could be continued up to the full system 
level. This will be done formally here. In order to make the method clearer, 
using our existing notation let us first evaluate the full set of equations for 
the undirected line graph with three nodes which we refer to as the open 
triple, depicted in Figure 2. Here we shall assume that the transmission rate 



parameter is r across both links and that the removal rate is 7. Firstly we 
write all of the single node equations for this network. From (1.1): 




2 



Fig 2. Open triple graph. 



(2.1) 



ih) 
ih) 



T{Sil2)-7{h), 
r{hS2)+T{S2h)-j{l2) 

r{l2S3)-l{h), 



and 



(2.2) 



(Si) 
{S2) 
{S3) 



T{hS2) - t(52/3) 

rihSs). 
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We also need to specify the following equations for pairs: 

(Sih) = T{SiS2h)-r{Sil2)-7{Sil2), 
(/1S2) = -T{IiS2h)-r{IiS2)-^{IiS2), 

{S2h) = -T{hS2h)-T{S2h)-^{S2h), 

(2.3) (/253) = T{hS2S3)-T{l2S3)-7{l2S3). 

Finally, at the triple level we have from the master equation (since the system 
has only three nodes): 

{SiS2h) = -T{SiS2h)-7{SiS2h), 

{I1S2I3) = -2T{hS2h)-2^{hS2h), 

(2.4) {118283) = -T {118283) -7(118283). 

In order to formulate the full system for an arbitrary graph, we introduce 
notation for the subsystem states. 

2.1. Notation for system and subsystem states. In general, our stochastic 
system (which we denote by T) comprises of P individuals, each of which may 
be in any of the 8, I or R states at any given time. In total, this corresponds 
to 3^ possible states. Denoting these system states by F", a G {1, 2, 3''^}, 
the probabilities for each state are given by the master equation: 

3^ 3^ 

(2.5) (F") = ^ (r^) - I] c^^" (r") , 

13=1 ^=1 

where a denotes a constant matrix of Poisson rate parameters. The master 
equation completely describes our stochastic system using a set of 3^ or- 
dinary differential equations. Our overall objective is to show that (1.3) is 
implied by the master equation when T represents a network with no loops. 

It is useful for us to define a general subsystem ipw comprising of r nodes 
in F indexed by vector W of length r: W = (W^i, W2, ...,Wi, ...,Wr), Wi G 
{1,2,..., P} where we can assume Wi < W2 < ... < Wj-. We assume that the 
network connections of the nodes of ipw Eire a subset of the connections of 
F. 

Let denote the state of subsystem ij^w where A = {Ai, A2, Ar) 
and Ai € {8,I,R} Vi € {l,2,...,r} is a sequence of 8,1 and R symbols 
of length r such that the state of node Wi is Ai. In terms of the notation 
of the previous section for subsystems of single nodes and pairs of nodes, 
we have 8i = ipf , li = ipj , 8ilj = tpfj, 8i8jlk = i^fji etc. We shall use 
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these two notations interchangeably. We shall also sometimes treat indexing 
vectors such as W as sets such that n £ W means that the node n is in the 
subsystem ipw- 

In general, although we can specify the states of each node with this type 
of notation, an important ambiguity remains because information about the 
network structure is not included. To remove this ambiguity, the notation 
should normally be used in the context of a sketch of the relevant network 
structure or where the network structure is clear from the context of its use 
(as in (1.1)). 

Let us now show how the differential equations of the different subsystem 
states can be formulated in general. 

2.2. Differential equations for subsystems. Here we obtain differential 
equations describing the rate of change of the state of any subsystem. Firstly 
we define two useful operators which act on subsystems. 

Definition 2.1. For the subsystem il^^r, if node Wk is either susceptible 
or infectious then: 



Remark. This operator changes the state of node in subsystem ipw 
to S if it is infectious. If node Wk is susceptible or removed then it leaves 
the state unchanged. 

Definition 2.2. A neighbour of node i is a node with a network link 
directed towards i. 

Definition 2.3. Ni denotes the set of neighbours of node i. That is: 



Definition 2.4. For the subsystem state V'ly of r nodes, a subsystem 
of r + 1 nodes can be generated as follows: Take k € {1,2, ...,r} and take 
a neighbour n of Wk outside of the subsystem with a network link towards 
Wk, i-e. let n € ^w^; n ^W. If Ak = S, then the generated (r + 1) -state is 
given by the generating rule: 




Ai...Ak-iSAk+i...Ar 

'w 



If node Wk is removed then: 




, Ai...Afc_ii?Afc+i...^r 



Tijj^oyjeN,. 



9Wk .,Wr,' 
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i.e. the subsystem is extended by an infected at node n which is connected 
towards W^. // Afc = /, then the generated subsystem state is given by: 



i.e. the subsystem is extended by an infected at node n which is connected 
towards and the state of node Wk is changed from I to S. 

To complete the definition, if = R then the operator g^y leaves the 
subsystem unchanged. We also assume that for any state where there 
is no link from node n to node Wk in the transmission matrix T, then the 
subsystem is also left unchanged. 

Remark. The generated order r + 1 system is obtained by replacing a 
susceptible or infectious node Wk in the original subsystem by an SI edge 
such that the S node of the edge is put in the place of the node Wk and 
where the / node of the edge is external to the subsystem. 

Definition 2.5. For any subsystem tpw of r nodes in state and 
a single node k G {1,2, ...,r} in state '4^^^' ^^^'^ define D^"" to have 
value 1 if the states of the node Wk match (Ak = a) and to have value zero 
otherwise: 



Theorem 2.1. The rate of change of the probability of a subsystem state 



The proof of this theorem is a rather long diversion and can be found in 
Appendix A. 

As an example of applying the theorem, we can use it to obtain the set 
of subsystem equations (1.1) by considering each equation in turn: 
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• If the subsystem is a single susceptible individual ipf , then r = 1 so 
k can only take the value k = 1 where Wi = i and Ai = S, reducing 
(2.6) to: 

(Vf) = - Yl TU^li)- 

The first term on the second line of (2.6) is zero because Tu = 0, and 
the last term is zero because Df^ = 0. 

• For an infectious individual ^pj we obtain: 

{4,0= ^2 TUi^^i)-l^{^l)^ 

where the first term on the second line of (2.6) is zero because Ta = 0. 

• If the subsystem is the pair ipfj then the sum over k is over k = 1 and 
k = 2 and Wi = i, W2 = j, Ai = S, A2 = I so: 

n=l,n^{i,j} 
P 

n=l,n^{i,j} 

where the first line corresponds to /c = 1 and the second to A; = 2. 

• If the subsystem is the pair then the sum is over k = 1 and k = 2 
where Wi = i, W2 = j, ^1 = 6' and A2 = S so: 

(^i;/) = - E T,^(^^)- E r.n(0' 

n=l,n^{i,j} n=l,n^{i,j} 

where both terms come from the first line of (2.6). 
We have therefore obtained (1.1) in a slightly different notation. 

2.3. The state space for a tree-graph. Here we build up a state space 
which is sufficient to describe a tree-like graph. We first make some defini- 
tions. 

Definition 2.6. An r -motif is a subsystem ofT comprising of r nodes 
and network links such that it forms a weakly connected network. 

Definition 2.7. An r-state is the state of an r-motif. 
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The state space that we need to consider is built up inductively from 
the states of single nodes by considering the infection process. Starting with 
the infected states of the single nodes i G {1, 2, P}, (2.6) shows that 
they depend on the 2-states i^fj, j G Ni as described by the generating rule 
(Definition 2.4). 

The differential equations for {t/jfj) in turn contain the 3-states ipfji: 
k G Nj and ipl^/j, k £ N^. The differential equations for the 3-states contain 
4-states and typically, the differential equations for r-states contain (r -|- 1)- 
states for r € {1, 2, (P — 1)}. This state generation process can continue 
until we reach P-states which can only depend on other P-states. 

Notice that this process always forms subsystems which are motifs and 
that the motif states can never include removed nodes. 

Proposition 2.1. For a tree-like graph, if the neighbours of the I nodes 
are all S in an r -motif, then this is true for all r + 1-motifs generated from 
this r-motif. 

Proof. This follows easily from the definition of the generating rule (Def- 
inition 2.4). □ 

Definition 2.8. Consider a tree-like graph and take the 1-motifs with 
I nodes: ip(,i G {1,2, ...,P}. The "basic state space" M is formed by these 
1-states together with the set of motif states that can be iteratively generated 
from them using the generating rule (Definition 2.4)- 

Remark. Due to the method of its construction, the state space M gives 
a self-contained system of differential equations, i.e. the time derivatives 
of the probabilities of each motif state can be expressed in terms of the 
probabilities of other motif states in the state space. An example in the case 
of the open triple is given by the motifs in (2.1), (2.3) and (2.4). 

Definition 2.9. Consider a tree-like graph and the 1-states: {ifj(,\(j^ : 
i G {1,...,P}} and the 2-states with SS, i.e. {V'f/ : i G {l,...,P},j G iVj. 
The "extended state space" M comprises of these motifs states together with 
the set of motif states that can be generated from them by repeated iteration 
of the generating rule. 

Remark. This extended state space is required to form the relevant clo- 
sure relations. Due to the method of its construction, it is also self-contained. 

Lemma 2.1. Let il)w ^ Then the neighbour of an I node is an S 
node in 
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Proof. Follows from Proposition 2.1. 



□ 



Lemma 2.2. For a tree-like graph, the equation for the probability of an 
r-state ip^r G M is given by: 



k=l 



n=l,niW 



(2.7) 



1=1 



Proof. For these states we have 1 - D^^ = 1 V/c € {1, 2, r}. Addi- 
tionally, when D^^ = 1, the first term on the second line of (2.6) never arises 
because Df^^ = 1 implies that an / is connected to an / node in r-state V'fy 
which contradicts Lemma 2.1. Therefore (2.6) reduces to (2.7). □ 

Let us now formulate the exact closure relations and prove our main 
result. 

3. Closure relation and proof of the main result. The exactness 
of (1.3) is straightforward to see provided that outbreaks of epidemics are 
always initiated with a single infected individual. We prove this first before 
considering the general case. 

3.1. Proof for single initial infected. When infection is initiated on a 
tree-like graph at a single individual, infection must always proceed in linear 
chains. Consequently there is no possibility of the state IkSjIi illustrated in 
Figure 3 arising because an infection initiated at either k or i must pass 




Si 

Fig 3. Shows a state which cannot arise on a tree-like graph where there is only one 
initially infectious node. 

through j to get to the other node. Furthermore, 

{Sjh) = {S^Sjh) + {liSjh) + {R^Sjh), 
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but since (liSjIk) = and consequently (RiSjIk) = 0, we have: 



reducing (1.1) to the following closed system: 



p 



{S^) 



P 



p 



Tjk{SjIk) - Tij{SiIj) - ^j{SiIj) 



P P 



{S,Sj) 



TikihSi) — ^ Tjk{SjIk). 



Similar arguments show that this can be written in the form of (1.3). 

More generally, this argument also applies to any tree-like graph where 
there is at most one network path by which any susceptible individual in 
the network can become infectious from the initial configuration of infected 
individuals. 

For a more general starting condition, we cannot use such a simple ar- 
gument. Before discussing the general proof for any tree- like network, we 
consider two very simple example networks which will serve to motivate and 
illustrate the method of proof. 

3.2. Proof for an open triple. Here we consider the case for the open 
triple depicted in Figure 2. The equations for the probabilities of the basic 
state space M are given in (2.1), (2.3) and (2.4). To form the relevant closure 
relations, we require the equations for the extended state space M formed 
by the equations for M together with (2.2) and: 



(3.1) {S;S2) = -T{SiS2l3) and {8283) = -t {1^8283). 



Our objective is to close the system at the level of pairs using (1.2), elimi- 
nating the need for (2.4), and show that the system remains exact. We note 
that the exactness of (1.3) can be proved in this case along the lines of the 
previous argument by considering each possible initial condition separately, 
however the approach discussed here will be more useful for understanding 
the general case. 
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We need to consider closures for the triples {I1S2I3), {I1S2S3) and {SiS2l3)- 
Let us consider the closure: 

Hlb2h) ~ TTTT • 

This is exact if a{t) = where: 

a{t) = {S2){IlS2l3) - {IlS2){S2h). 

Taking the derivative of a with respect to time gives: 

ait) = {S2){hS2l3) + {S2){IlS2l3) - {hS2){S2l3) " {hS2) {S2I3) ■ 

Substituting the relevant derivatives in from (2.2)-(2.4) and cancelling terms 
reduces this to: 

a{t) = -2(r + 7)a(t), 

so: 

a(t) = a(0)e-2("+^)*. 

Now it is easily verified that provided the system is initiated in a specific 
system state then a(0) =0 (the general proof of this statement in provided 
by Lemma 3.2 in Section 3.5). Consequently a{t) = for all t > and the 
closure is exact. 

By symmetry, it will suffice to consider one of the remaining two triples 
in (2.3). We wish to show that a{t) = where: 

(3.2) a{t) = (52) (5i 52/3) - {SiS2){S2l3)- 

Here it is necessary to also use (3.1) for pairs of type SS in the extended state 
space. Repeating exactly the same analysis as above shows that a{t) = for 
this closure as well, demonstrating that (1.3) is exact for this system when 
the system is initiated in a specific system state. 

3.3. Proof for a star graph. We now consider the case of the undirected 
star graph with P = A shown in Figure 4, where again we assume that the 
strength is the same across each network link and is denoted by r and the 
removal rate for each node is 7. Writing down the equations of the extended 
state space, there are two types of closure which need to be proved: one for 
the 5 — 5 — / triples and one for the I—S—I triples (see (1.1)). The graph has 
three triples ((1, 4, 3), (2, 4, 3), (1, 4, 2)) , but it is sufficient to prove exactness 
for one of them. Hence we want to prove the following two relations: 

(S4) (51/354) = (5i 54) (1354), 

(3.3) (54)(/i/354) = (/i 54) (1354). 
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Fig 4. Star graph with P = 4 nodes. 
For brevity, we adopt the alternative notation: 

(V'f)(^f3^4)-(V'f,f)(V'li) = 0, 
(V'f)(ff£4)-«i>(C4) = 0. 

We introduce: 

«i = (V'f)(V'fr4)-(V'i^f)(V'li)- 

By differentiating this, substituting in from the process equations and group- 
ing terms we obtain: 

(3.4) di = — (r + 7)ai — ra2 — tos — ra^, 

where: 



02 








= (^2^^4>(^lT4) 




a4 




-(^iT4)(^li) 



Differentiating 02 we get: 

(3.5) 02 = -2(t + 7)0:2 - TOs - rag, 
where 

«5 = (^f,^^)(^lT4)- (^1^4X3^), 
«6 = (V'(^4)(^fi!3"4)-(^f|)«2!f,4). 

The derivatives of 03 and 04 can be obtained similarly. 
Differentiating 05 we get: 

(3.6) ds = -3(r + 7)05 - ra7 - rasj 
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where: 




The derivative for ag can also be obtained. Finally, differentiating a-j and 
as we obtain: 



To conclude the proof of the exactness of the closure, we first assume that the 
initial state is not mixed; that is, one of the 3^ = 81 possible configurations 
has probability 1 at t = 0. Then it is easy to see that aj(0) = for all j 
(see Lemma 3.2 in Section 3.5 for a general proof). Hence the differential 
equations for a-j and as show that a'j(t) = and a^{i) = for all t > 0. The 
differential equation for as then implies that as = Vi > (and similarly for 
ae). This implies that a2 = (and similarly as = a4 = 0). The differential 
equation for a\ shows that ai = which is what we wanted to show. The 
other triple closure in (3.3) can be proved similarly. 

Remark. In fact, we have proved several closure relations aia2, ...,a8. 

The closure relations each consist of two pairs which are visualised in 
Figure 5. For reference, we refer to these as the left pair and the right pair 
referring to their position in this figure. Looking at these closure relations, 
we can form two observations: 

1. For a given node i, the number of times it appears as Si is the same in 
the left and the right pair, and similarly with the number of times it 
appears as Ij. For example, with as, node 1 (the left node) has one / 
and one S for both pairs and node 4 (central node) has two 5's in both 
pairs. For as, the number of I's at nodes 1,2,3 and 4 is (1, 1, 1,0) in 
both pairs and the number of 5's is given by (1, 0, 0, 2) in both pairs. 

2. Any SI pairing on the left appears exactly the same number of times 
on the right. For example in ay, I\S/^ appears once on the left and 
once on the right and 1-^,3 appears twice on the left and twice on the 
right. Observing that only SS pairs and IS pairs appear in the closure 
relations, a consequence is that SS pairs also have this property. 

These observations will be of key importance for developing the general proof 
in the following two sections. 



(3.7) 



dy = — 4(r + 7)a7 



and 



ds = - 4(t + 7)as. 
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Fig 5. Each box illustrates the relevant node states for the four parts of the closure rela- 
tion in the equation above it. The node states on the left and the right correspond to the 
two terms m the closure relation. The node numbers correspond to the same positions as 
Figure 4- 
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3.4. General closure relations. In general, to show that the closure re- 
lationship (1.2) is exact for the tree graph, we need to show that a = 
where 

a = {Bj){A,BjCk) - {AiB,){B,Ck). 

Our proof of this is via induction using a sequence of closures analogous to 
the proof in the case of the star graph in Section 3.3. 

We shall consider many closure relations. In general we specify that they 
are composed of two pairs of motif states {ipw-''^x) ^^"^ i.'^^Y^'^z) ^^"^ thai 
the closure is exact if a = where 

« = (V'i^)(V'l)-(V'?)(^f). 

We formalise the observations we made about the closure relations for the 
star graph at the end of Section 3.3 by defining what we term "compatible 
pairs" . 

Definition 3.1. ipf c V'h' ^3 — ^ ^"-^ = 

da C ^ 3ii, j2 s.t. Wj, = i, Wj^ = 12 and Aj^ = ai,Aj^ = 03. 

Definition 3.2. Two pairs of motif states (V'w^iV'x) '^'^^ (V'yjV'z) ^'^s 
called compatible pairs if the following conditions are met: 

• CP(i) i/j^ C Viy or Tpf c V"! ^ C C ^pY or V'f C V'f • 

• CP (a) C Viy and Cipx ^ ^"■^ '^i V'f • 

• CP (Hi) C iJ^ or C V'l ^ C or C • 

• CP(iv) C i;^ and i^lf^^^ C V'l ^ C and C • 

• CP(v) Same as CP(iii) and CP(iv) hut with SS pairs. 

Definition 3.3. Let be an r-state and tp^ a q-state. Then the 
order of the pair (ipw^'^x) defined as r + q. 

Proposition 3.1. If{7p^,ipx) and (ipyjipz) are compatible pairs, then 
their order is equal. 

Proof. Follows from CP(i) and CP(ii). □ 

Proposition 3.2. For a tree graph, applying the transformation hi to 
each of the four motif states in compatible pairs that contain node i generates 
compatible pairs. 

Proof. The transformation satisfies CP(i) and CP(ii) because it replaces 
Tpf = tpl with Tpf = Tpf which does not alter the form of the conditions. 
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The transformation satisfies CP(iii) and CP(iv) because all IS pairs where 
i is the infected individual are removed by this transformation. New IS 
pairs cannot be created by the transformation since this would require II 
pairs which are prohibited for tree graphs by Lemma 2.1. CP(v) is satisfied 
because the transformation leaves existing SS pairs unchanged and created 
SS pairs result from existing IS pairs so are balanced on each side. □ 

3.5. Proof of the main result. 

Lemma 3.1. Let (V'w'; V'x) ^'^^ (V'yi V'z) compatible pairs or order R 
and iiy^,'4)x,^Y^'^z ^ -^et- 

«o = (V'^)(V'l)-(^?)(^f). 

Then: 

m 

(3.8) «o = X! "^P^P + "^oao 

p=i 

where each ap can be expressed as: 

where ip^^ip^ and o'^s compatible pairs of order R + 1 and where 

Co, Cp are constants and m is an integer denoting the number of terms in 
the summation. 

Remark. This is a general statement of the forms of (3.4)-(3.7) in the 
star graph example. 

Proof. Take the derivative of oq: 

(3.9) do = + - ) - )• 

We consider the terms associated with removal, transmission terms of order 
R and transmission terms of order R + 1 separately. Firstly, from (2.7), this 
derivative contains the following terms associated with the removal process: 

ki k2 

+ E Dg'^n^ ^^y) (V'f ) + E (V^?) (V^f ) 
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where the sums over /ci , /c2 , /ca , /c4 are over ah nodes in the motifs ipw , ipx , V'y i V'z 
respectively and where 

V = E K'^w,. + E = E <7y.3 + E ^f/7z., 

ki k2 fcs fc4 

is easily seen to follow from CP(i) and CP(ii). 

The right hand side of (3.9) also contains the following transmission terms 
with motifs of order R: 

ki h 

-E<'Erx..^..Af(V'i^)(^l) 

+E^^3"E%3^.A?(^?)(^^) 
+E<"Er..,z,,A?^(^?)(^f) 

= -wao, 

where: 

- = E<'ErH..,H.,Af + E<"E^^..^..Af 

= E Dg' E Tn^y.M' + E <^ T.Tz..z^dFJ 

fcs '3 ki I4 

and where the sums over /ci, A;2, /cs, A;4, /i, ^3) ^4 are over all nodes in each of 
the relevant motifs. This follows from CP(iii) and CP(iv). Hence the removal 
terms and transmission terms of order R contribute coao to the derivative 
of ao where cq = —v — w. 

For transmission terms with motifs of order i? + 1, consider the term 
{ipw)i'^x) ™ (3.9). This gives rise to the following terms in the derivative 
of ao'. 

r~\AS ^ 
ki n=l,n0V 

To prove the lemma, it is sufficient to show that each term in this sum can 
be paired uniquely with a term in {iPy){'^z) '^^ ™ (V'y')(V'z)' such that the 
difference of these terms forms ap. By symmetry this is a one-to-one pairing 
establishing that each term of order i? + 1 on the right hand side of (3.9) is 
accounted for exactly once in the sum of Op. 

imsart-aap ver. 2012/08/31 file: Exactequations.tex date: December 11, 2012 



EXACT EQUATIONS 



21 



Let us take an element from the sum by choosing a node W^, w € 
{l,2,...,r} and an outside neighbour node n € Nyy,^. This neighbour can 
either be in ipx or outside. 
Case 1: = S 

Consider first the case where is a susceptible node, resulting in the 
following term in the sum: 

where we can identify Cp = —T^^n- 

We have Ay^ = 5", ?i G ^^Ww ^^'^ ^ Let us now identify a term in 
{'4>y){'4'z) + ('V^y)('^z) *o form compatible pairs. According to CP(i) and 
CP(ii) we can assume without loss of generality that Wyj € Y , i.e. 3j/ : = 
Wyj and Cy = S (see Figure 6). There are two subcases: 
Subcase 1.1: n ^ Y 

When n ^ Y, either n G Z or n ^ Z and these are shown by solid and 
dashed lines respectively in Figure 6a. By CP(i), n Z 4=^ n ^ X since 
n ^ so the solid lines match on the left and right pairs as do the dashed 
lines. The corresponding term must therefore be —TYyn{9Yy(.'4^Y))i'^z) ~ 
'^pidYyi'^i^v)) ^"^z) irrespective of whether n G Z or not. Hence: 

where {9w^('^w)^'^x) {9Yyi'^Y)''^z) easily seen to satisfy the def- 
inition of compatible pairs since the extra node is n which is / in both 
pairs. 

Subcase 1.2: n £Y 

If n G y, then the edge n ^ Yy is an SS or IS edge in C. By CP(iii), 
CP(iv) and CP(v), it is also the same edge in B because n ^ W. Hence 
3x : Xx = Wu, and Bj. = S. By CP(ii), £ Z is also true where 3z : 
Zz = Ww and = S. This is illustrated in Figure 6b. We therefore have 
the corresponding term -7z,n('0y ) (^z, (V'f )) = Cp(?/'y)(ffz,(V'f )) and: 

«p = - (V'?>(ff^.(^f )), 

where again, the relevant pairs are seen to satisfy compatibility. 
Case 2: Ayj = I 

So far we have proved the existence of oip when A^ = S, n £ Nw^ , n ^ W. 
Now we have to show Op can be defined when A^j = I, n £ A^iy,„, n- ^ W. 
In this case, the motif generating rule firstly changes A^ = I to A^ = S 
and then applies the same generating rule as if A^ = S initially. From 
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a) Subcase 1.1 




b) Subcase 1.2 




Fig 6. Each circle refers to one of the motif states -ifj^ ,^1}x,'4^y ^ "4^ z specified to the top left. 
The position of the relevant node states with respect to the motif states are then illustrated, 
a) Subcase 1.1 ( n^Y). b) Subcase 1.2 (n e Y). 
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Proposition 3.2, applying the transformation to compatible pairs of order R 
produces compatible pairs of order R in the case of the tree graph. After 
this transformation, the argument runs identically to case 1. This completes 
the proof of Lemma 3.1. □ 

Lemma 3.2. Assume that the initial condition is not mixed, i.e. 3A € 
{I,S}^ s.t. {ipi2 p) = 1- If {i'w^'^x) '^^'^ (V'y-j'^z) ^'''^ compatible pairs, 
then for any graph, - (V'y)(^f ) = t = 0. 

Proof. Assume that = 1- Then by CP(i) and CP(ii), for all 

Tpl C 1p^r^ we must have ijjf C and/or ipl C V'z- This is also true for 
all ipi C ipx ^^"^i by the symmetry between the compatible pairs, it follows 
that {ipyYi'^z) = 1- Similarly, it follows that {i>y/){'^x) = implies that 

(V'?)(^f) = o. □ 

Theorem 3.1. Let us assume the following: 

• The graph is a tree (no loops) 

• The initial condition is not mixed. 

• (V'w/jV'x) {'^Yi'^z) ^'''^ compatible pairs. 
Then {tp^){tp^) - (^y)(^f ) = for all time t > 0. 

Proof. We prove the theorem by induction according to the order of the 
closure. This is analogous to the proof for the star graph in Section 3.3. 
Step 1 

If the closure is of order 2P, then it is exact. More precisely, if '^l>^r, il'x^ 
ipy and tp^ are P-states, then (3.8) does not contain the summation terms 
and becomes: 

do = coao- 

Since we start from an initial condition that is not mixed, we have (by 
Lemma 3.2) ao(0) = ^ ao{t) = Vt > 0. 
Step 2 

Assume that the theorem is proved for compatible pairs of order R+1. We 
prove that it is true for compatible pairs of order R. Applying Lemma 3.1, 
we have: 

m 

do = X! + "^oao- 

p=i 

According to the induction condition, Op = Vp because these are com- 
patible pairs of order R+1. Therefore do = coaQ. From Lemma 3.2, ao(0) = 
so ao(*) = Vt > 0. Since we have proved the result for compatible pairs of 
order 2P then we have completed the proof of the theorem. □ 
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The lowest order compatible pairs are of order four. The closure rela- 
tion corresponding to these pairs is formulated in the following important 
corollary. 

Corollary 3.1. Under the assumptions on the graph and the initial 
conditions in Theorem 3.1, we have the special cases: 

(^f)(^f;f.) = (^i:f)(f.) 

for all i S {1, 2, P} and for all j S Ni, k ^ Nj, i ^ k; 

for all i € {1, 2, P} and for all j, k S Ni, j ^ k. 

This corollary is Theorem 1.1 expressed in a different notation. 

4. Application to some graphs which are not tree-like. To com- 
plete this work, we make a final observation which shows that the pair-based 
model can sometimes provide an exact representation of the infectious dy- 
namics of graphs which are not strictly tree-like. We first make two defini- 
tions which can be understood with reference to the examples in Figure 7. 

Definition 4.1. A reduced representation is a graph which is constructed 
from the initial transmission network and the given initial conditions by re- 
moving transmission routes which cannot carry infection dynamics. 

Definition 4.2. An independent segment is a region of a graph that is 
only connected to other regions via nodes which are initially infectious. 

Theorem 4.1. Given SIR dynamics on a transmission network with in- 
fection and removal governed by Poisson processes and given an unmixed 
initial state of the system, if every independent segment of the reduced rep- 
resentation is tree-like, then applying (1.3) to this representation exactly 
generates the expected infection dynamics on the original transmission net- 
work. 

Proof. By definition, the infection dynamics of the system remain un- 
changed after the removal of edges which cannot support infection dynamics. 
Additionally, the infection dynamics of any independent segment are inde- 
pendent of the dynamics on the rest of the graph because there is no process 
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Original grapti Representation grapli 




I I 




Fig 7. The graphs on the left are the initial transmission networks where the initially 
infected nodes are indicated by the symbol I. The graphs on the right-hand side are the 
reduced representation graphs where the cuts for independent segments which occur for 
cases b and d are indicated with dashed lines. The tree-like structure of the graphs on the 
right-hand side shows that applying the pair-based model to these graphs generates an exact 
representation of the infection dynamics on the original system. 
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that allows influence across the initially infectious nodes. If the resulting rep- 
resentation graph is a set of tree-like segments, then since (1.3) is an exact 
representation of the dynamics on each independent segment, solving (1.3) 
on the reduced representation graph is equivalent to the infection dynamics 
on the original transmission network. □ 

Figure 7 shows some graphs and the associated representation graphs 
where the dashed lines indicate the boundaries that separate independent 
segments. For each of these examples, the solution of (1.3) on the repre- 
sentation graph exactly reproduces the expected infection dynamics of the 
original system. 

This suggests that the accuracy of the pair-based model could be increased 
by first generating the representation graph for the particular network and 
initial conditions prior to numerically solving the pair-based model. 

5. Discussion. We considered the pair-based variant of the subsystem 
approach to constructing epidemic models on networks [17, 18]. We proved 
that for SIR dynamics on fixed tree-like networks with exponentially dis- 
tributed transmission and removal processes, the pair-based model provides 
an exact determination of the infection probability time course for each indi- 
vidual in the network. We also showed that the dynamics of some networks 
with clustering can also be represented exactly by the pair-based model 
under specific initial conditions. 

This represents the first viable deterministic model of epidemic dynam- 
ics on finite heterogeneous systems which exactly matches the underlying 
stochastic processes. We use the qualifying term "viable" here to exclude 
models based directly on the master equation because these are too numeri- 
cally intensive to be practical for systems of reasonable size and complexity. 
Additionally, "heterogeneous" excludes systems with significant symmetry 
which may be employed to obtain exact representations in very specialised 
circumstances [10, 19]. 

We note that effective degree models can generate very good agreement 
with stochastic simulation [3, 12] as do the PGF or edge-based compartmen- 
tal modelling methods [13, 14, 22], although exact correspondence has not 
been proven. For some idealised networks, including fully connected net- 
works and some configuration networks [22], convergence to the expected 
value can be shown in the infinite population limit [3, 5]. However, these 
models have a large measure of homogeneity, and convergence only occurs 
for infinite populations. 

It is intuitively understood that clustering is at the root of problems with 
models based around closures at the level of pairs [9]. Previous analysis [18] 
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attributed the failure to anomalous terms which emerge in subsystem equa- 
tions when differentiating closure approximations based around the statisti- 
cal independence of individuals. Here, repeating similar analysis for a closure 
at the order of pairs (1.2) in the context of tree- like graphs, these anomalies 
do not arise and we are able to prove that the closure is exact via induction. 

In principle, models based around subsystems at the order of three nodes 
or higher could be constructed. The next higher-order model would require 
obtaining a closure which is able to preserve correlations between triples, 
and similarly for higher orders. This leads to an interesting theoretical ques- 
tion for future analysis: does the hierarchy of exact order-by-order models 
suggested by Sharkey in [18] exist, and if so, what form should the closure 
approximations take at each level? We conjecture that exact closures of a 
similar nature to those considered here are possible for networks with more 
structure, given that the order at which the closure is performed is guided 
by the network structure. For example a network with loops of at most size 
three, and with no additional structure should be closed so that triples are 
conserved, or at the levels of quadruples. This seemingly heuristic obser- 
vation follows from the compatible pair argument presented in this paper. 
Furthermore, the problems of identifying the smallest motif where compati- 
ble pairs can be found is directly related to a graph colouring problem which 
we will expand upon in future work. 



The proof of Theorem 2.1 is analogous to the proof of the single and pair 
equations by Sharkey (2011) in [18]. In what follows, summations over Greek 
indices a,/3 are assumed to be over all 3''^ possible system states. First we 
make some definitions. 

Definition A.l. For a system T in state a and a single node iofTin 
state a we define: 



denoting whether or not the specified single node state matches the system 
state. Note that this is just Definition 2.5 applied to the full system. 

Definition A. 2. 



APPENDIX A 




otherwise 





1 if the states of all individuals in V are the same 

for and for F^ except for ipj which may change 
otherwise 
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Proposition A.l. For all a,i: 

a 

where the summation is over all possible states available to node i. 



Proof. Statement that for a given system state P", or subsystem state 



each node must be in a unique state. □ 



Proposition A. 2. For all f3,i,a: 



Proof. Statement that there is only one system state which is identical 
to r'^ except that node i is in state V'f • D 

Proposition A. 3. For any subsystem ip^ and Mk G {1, 2, r}; 
for all a, a. 



Proof. Proposition is true when Dy^ = 0. When Dy^ = 1 we have: 



□ 



Proposition A. 4. For any subsystem ^1;^^ andMk G {1, 2, r}; 

r r 
^wSWk 11 ^W^ ^Wj - 11 

for all P, a. 

Proof. Proposition is true when: 

r 



n <7=o- 

j=^J¥=k 
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From Proposition A. 2 there must be a single state for which D^^(^^ = 1, 
otherwise it is zero. When: 



n D^'=^^ 



we must also have (for the state when Dy^^^^^ = 1): 



n <^ = i 



because only site ipw^ can change state during this transition, establishing 
the proposition. □ 

We can now use these propositions to prove Theorem 2.1: 

Proof. We have that: 



a 1=1 



Taking the derivative of this with respect to time and substituting in the 
system master equation (2.5) gives: 



a i=l 



(A.l) 

From Proposition A.l: 



i=l ap i=l 



1 



ED 

ar 



aar 

Wr 



bi 



'Wr 



br 



= EEn^^^^^ _ 

fe=l a^bk j=l 

Multiplying the right of (A.l) by this gives: 



5:."^(F^)n<^EEn^^^''^' 

0/3 i=l fc=la(;6fej=l 



a/S i=l k=l a^bk i=l 
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This can be simplified using the fact that cj"'^ = whenever the state of 
the subsystem ipw differs by more than a single individual ipWkJ ^ S {l...r} 
between states and which means that aj = bj = Aj for j ^ k: 

a/3 i=l k=lj,j^k afcfefe 

-E-'"(nfK^E n <^<tE^^:< 

a/S i=l k=lj,jjtk afcfefe 

= E-"''(r'^)E<^ n <^<t E^s?;<: 

a/3 k=l o.kbk 

r r 

a/3 fc=l jjVfc afcbfe 

where the last equality follows from D"^,^ D'^,^ = . 

For SIR dynamics, we can do the summations over and hk'- 

{i^^) = E-"'(r'^)E<^- n <^<t« 

a/3 k=l j=l,jjtk 

a/3 A.-=l j=ljyfc 

(A.2) -E^'"(r")E<^ n <^<t«- 

a/3 fc=l i=i,jyfe 

Now we introduce the relevant terms in the transition matrix at the level 
of the system: 



r _P r 

r P r 

r r 

where these equations are designed so that they are satisfied for any combi- 
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nation of a, /3, k. Substituting these into (A. 2) gives: 

r r P r 



^3 

(xfi k=l j=l^j^kn=l j=ljz/zk 

-E(r")E<^ n E^H^.-«^^?'c^: n <^<t 

a/? k=l i=l,i^fcn=l j=lj^k 

r 



-E(r")E<^ n 7H^.«cs^ n <^<^ 

a/? fc=l 1=1, i^A; j=^d7^k 

We can rearrange the summation order: 

(^i^) = EE^H^."E(r^)<^fE<^<C^ n <^<t 
-EE^H^.-E(r")<^<^n^E<c^: n <^<t 

fc=l ji=l o /3 j=ijjkk 

-E7iy.E(r")<^<E<c^: ri <^<t' 

fc=l a IS j=l,jf^k 

and apply Proposition A. 3: 

(K) = E^.'''E^H..nE(r^)<^n''E<C! n <^<t 

fc=l 71=1 /3 a j=ljjLk 

-E^.^"E^^.nE(n<i^n"^E<c^. fl <^<t 

k=l n=l a /3 j=ljjLk 



E^.^V.E(r"XE<C. n <^<7- 

fe=l a 13 j=l,j7^k 



Applying Proposition A. 4 gives: 

fc=l n=l /3 j=l,j^k 

fc=l n=l a j^ij^f^ 



k=l a j=ijj^k 
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Breaking up the sums over n in the first two hues depending on whether the 
node n is internal or external to the motif ipw gives: 



k=l n^W P j=l,j¥'k 

r r 

PA, 



k=l n£W P j = lJjLk 

k=l n^W a j=l,j7^k 

k=l nGW a j=l,j^fc 

k=l a j=l,3^k 



Lines one and three can be immediately recognised as the generating rule 
(Definition 2.4). For n^W and n € -/Viy^: 

and if n ^ then Tvi/j.n = 0. 

Line two requires that n € W . Let / G {1, 2, r} and Wi = n. Then: 

tTw.w^j:i^')DZD^j n 

1=1 p j=ijVfc 

= t Tw.w, E(r^x«t n <t 

= E^H..H.,A^^E(r^)« n <t' 

where the last equality follows from Proposition A. 3. Using the definition of 
hwhii^w)' ^^^^ becomes: 

r 

1=1 
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and similarly for line four. 
We obtain: 

k=l n0V 
k=l 1=1 
k=l n(^W 

- E E Df'Tw.w, {hw, (K)) 

k=l 1=1 
k=l 

where the operator on the fourth line is superfluous but allows us to 
write the equation in the form of (2.6). □ 
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